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Abstract 

An effective Hamiltonian for the ferroelectric transition in PbTiOs is constructed from first- 
principles density-functional-theory total-energy and linear-response calculations through 
the use of a localized, symmetrized basis set of "lattice Wannier functions." Preliminary 
results of Monte Carlo simulations for this system show a first-order cubic-tetragonal tran- 
sition at 660 K. The involvement of the Pb atom in the lattice instability and the coupling 
of local distortions to strain are found to be particularly important in producing the be- 
havior characteristic of the PbTiO^ transition. A tentative explanation for the presence 
of local distortions experimentally observed above T c is suggested. Further applications of 
this method to a variety of systems and structures are proposed for first-principles study 
of finite-temperature structural properties in individual materials. 
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electrics, February 1995. 



1. Introduction 



Structural phase transitions in perovskite structure compounds have long been of ex- 
perimental and theoretical interest This class of materials exhibits a wide variety of 
transition behavior and low temperature distorted structures, depending on the individual 
compound. From first-principles calculations, quantitative information on the structural 
energetics can be obtained to understand these chemical trends. Furthermore, this ap- 
proach yields a microscopic description of the origin and nature of lattice instabilities in 
these materials. 

Success in applying the first-principles approach to perovskite oxides has been achieved 
only rather recently, due to the presence of oxygen and transition metal elements, the num- 
ber of atoms in the unit cell, and the high accuracy required for a meaningful computation 
of the delicate lattice instabilities. The 1992 LAPW study [§ of BaTi0 3 and PbTi0 3 has 
been followed by investigations in a number of systems using various implementations of 
density-functional-theory (DFT) methods @] 0] @- In this work, we combine 
total energy and electronic bandstructure calculations with the DFT perturbation method 



||1C|| | fn| ||12| The latter provides a means for computing force constant matrices at 

arbitrary wavevectors q without the need for supercells, thus facilitating the exploration of 
lattice instabilities throughout the first Brillouin zone. Furthermore, valuable information 
about polarization density in these systems can be obtained from direct computation of 
the Born effective charge tensors and electronic dielectric constant e^. 

The highly demanding nature of these computations has the consequence that direct 
molecular dynamics or Monte Carlo simulations of temperature-driven structural transi- 
tions in perovskite oxides using DFT forces and total energies is completely impractical. 
Indeed, the investigations mentioned above do not attempt to map out a full configura- 
tional energy surface, but rather focus on those aspects felt to be most relevant to the 
description of the structural transitions: (i) lattice instabilities at quadratic order; (ii) 
anharmonicity at q = 0; and (iii) coupling of unstable modes to strain. We have developed 
a method based on the construction of a "lattice Wannier function" basis jnj which uses 
this first-principles information for an individual material efficiently to derive an effective 
Hamiltonian in the low-energy subspace of the ionic displacement space relevant to the 
temperature range of the transition. This effective Hamiltonian can be regarded as a real- 
istic generalization of the $ 4 models that are widely used in studies of generic ferroelectric 
transition behavior [|15] [|16| [T7[] A Monte Carlo simulation of this effective Hamiltonian, 



by construction, will quantitatively reproduce the T 7^ behavior of the particular system 
under consideration. 

To illustrate this approach, we present preliminary results of its application to the 
ferroelectric transition in PbTiO^. This compound has the cubic perovskite structure 
at high temperatures, undergoing a first-order transition at 763 K to a tetragonal phase 
produced by relative displacement of the cationic and oxygen sublattices resulting in a 
nonzero uniform polarization density. In the first-principles study of this transition, we 
obtain results for the transition temperature, the order of the transition and the nature 
of the paraelectric phase. The role of different contributions to the effective Hamiltonian, 
such as strain coupling and the q dependence of the lattice instability, can be examined. 
The significant differences from BaTiOs and SrTiOs |19| become evident. Although 
PbTiO^ has been described as a "textbook example of a displacive ferroelectric transition," 
H recent experimental evidence and our first-principles results suggest that the apparent 
simplicity of this behavior is quite nontrivial. 



2. Derivation of the Effective Hamiltonian 



The derivation of the effective Hamiltonian [14] is here summarized briefly. The start 



ing point of the analysis is the classical partition function 

Zoc J \du i )exp{-pHiat{{uj))) (2.1) 

where, as a result of the Born-Oppenheimer approximation, only the ionic degrees of 
freedom {uj} appear explicitly, the index j running over all ions in the system. The Taylor 
expansion of the full lattice Hamiltonian around a high-symmetry reference structure, 
such as the cubic perovskite structure, can be grouped into the harmonic terms TCh and 
higher-order terms in the displacements Tianh- Diagonalization of the harmonic part Tih 
is straightforward, and the results can be depicted using the dispersion relation of the 
type shown in fig. |l|. This suggests a natural decomposition of the full ionic displacement 
space into two subspaces: (i) Aq, composed of the branches of the dispersion relation that 
contain unstable modes characteristic of the structural transition, and (ii) A s , composed of 
the remaining branches, for which typically all modes have uj 2 (q) > 0. This decomposition 
permits the following approximation to the lattice Hamiltonian: 



Hiat « n h (A ) + n anh (A ) + n h (A s ) 
3 



(2.2) 



where only the anharmonic terms in the unstable subspace A , crucial for stabilizing the 
crystal, are included. As a result, the A s integral in the partition function becomes trivial, 
making no contribution to the temperature dependence of the average structure, and the 



An explicit form for the effective Hamiltonian 7i e //(Ao) is obtained by making a 
Taylor expansion in coordinates of the subspace Ao. The simplicity of this expression can 
be maximized by a proper choice of the corresponding basis. In analogy to the tight- 
binding representation of a subset of electronic eigenstates in a crystal, we choose a highly 
symmetric localized basis for A , where each basis vector Wi\ is associated with a Wyckoff 
position in the unit cell Ri and transforms according to an irreducible representation of 
its site symmetry group, the Wyckoff positions and irreps being determined by the space 
group irrep labels of A . These localized basis functions can be related to the Bloch vectors 
e^A, which diagonalize 7-^(A ), through an expression of the form Wi\ = J2^e iq ' Ri e^\ and, 
by analogy with the electronic case, are given the name "lattice Wannier functions." With 
this basis, an ionic configuration in the subspace Ao can be specified by the corresponding 
coordinates {^a}. 

3. Effective Hamiltonian for PbTiOs 

The first-principles method for the construction of the effective Hamiltonian for an 
individual material proceeds by the following steps. (I) The subspace A is identified, 
a symmetry analysis performed, and approximate Wannier basis vectors constructed by 
fitting to force-constant matrix eigenvectors at high-symmetry ^-points. (2) The form 
for the effective Hamiltonian is obtained by a Taylor expansion in symmetry-invariant 
combinations of the Wannier coordinates {^a}- The expansion is truncated to include a 
relatively small number of physically important terms. (3) Finally, the set of parameters 
is determined from first principles by fitting 7i. e ff to the results of selected total energy 
and linear response calculations, using the explicit correspondence, obtained from the 
construction in step (1), between the Wannier coordinate description {^a} and the actual 
ionic displacements {uj}- For each material, overdetermination of the parameters through 
additional independent first-principles calculations is recommended to establish the validity 
of the truncations made in step (2). 



partition function reduces to Z ps J A exp(— /?7i e //(A )) where 
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3.1. Construction of Wannier basis 

The relevant mode for the ferroelectric transition in PbTiO^, which freezes in to 
produce the low-temperature tetragonal structure, is the three-fold degenerate zone-center 
Ti5 mode. For a complete identification of the unstable subspace Ao, DFT linear response 
calculations were carried out at T, R, X and M (in the Brillouin zone of the simple cubic 
lattice) |HJ. The resulting low energy eigenmodes T^, i?i 5 , X' 5 , M' 2 and M§ require that 
the Wannier basis vectors Wi\, A = 1, 2, 3 are situated at the Wyckoff position occupied by 
the Pb atom and transform according to the vector irrep of the site symmetry group Oh- 
This results from the significant involvement of the Pb atom in the lattice instability 



and deviates from what would be expected if the transition could be viewed as resulting 



from the "rattling" Ti ion picture developed by Slater for BaTiO^ |p2[ , in which case the 
Wannier basis vectors would be situated at the Ti atom positions. 



Following the procedure described in detail in Ref. ||14||, the Wannier basis vectors 



are parametrized by identifying the displacement patterns centered at the Pb atoms which 
transform according to the vector irrep of Oh, starting from the innermost coordination 
shells. In the present construction, we include displacements of the central, first and second 
neighbor Pb ions and the nearest Ti and O shells, for a total of nine parameters (fig. |2|). 
The values obtained by fitting to the normalized eigenvectors calculated at Tis, i?is, M' 2 
and show a rapid decay with increasing distance from the central atom, as expected. 
Futhermore, the X§ eigenvector predicted from the parameters obtained from the fit at the 
other high-symmetry points is in excellent agreement with the normalized X§ eigenvector 



calculated from first principles. The results are given in detail in Ref. [[20 



3.2. Form of the effective Hamiltonian 

With this choice of basis, the system is described as a set of three-dimensional vectors 
at the sites of a simple cubic lattice, and 7i. e ff is expanded in symmetry-invariant 
combinations with respect to 0\, the space group of the cubic perovskite structure. In 
the present work, intercell interactions are included up to quadratic order only. For first, 
second and third neighbors, the most general quadratic interactions allowed by symmetry 
are included: 



J2 E MS ■ ^(5$ • ^) + Mii ■ £i(d) - (S ■ d)(&d) ■ d))] (3.i) 



+ E E ■ d){ii{d) ■ d) + b T i(£i ■ di)(S$ ■ di) + M& ■ 4)(&(d) • 4)] (3.2) 

1 d=nn2 

+ E E ■ ^(S(^ • 4) + c r(fi • - (S ■ <?)(£(<?) • d))], (3-3) 

1 d=nn3 

while beyond third neighbor we use a dipole-dipole form parametrized by the mode effective 
charge Z and the electronic dielectric constant e^: 

^ ^ (g ) 2 (g ■ - 3(g ■ <?)(g(<j) ■ d)) _ (3 4) 

i d - e °° Ml 3 

Terms in the onsite potential, depending only on values of £j at a single i, include isotropic 
terms up to eighth order in and full cubic anisotropy at fourth order: 

E( a i£i 2 + B &\ 4 + c (ti + 4 + A) + D &\ 6 + E \&\ 8 )- (3-5) 

i 

The energy associated with homogeneous strain, specified by the tensor e a p, (a, (3 = 
x,y,z), and its coupling to the local distortion £j, is included to lowest nontrivial order: 

+ y Cl1 1^ e ^ + ~2 12 e «« e /^ + "j * 4 2^ e «/3 + ^/ eaa 

+£/o(E ea °) E i^' 2 + 91 E( eaa E^«) + ^ 2 E 

a i a i ct<(5 i 

For our calculations in PbTiOs, we take e a p = at ao=3. 96883 A, the lattice constant 
measured for the cubic phase just above T c . We have further generalized this expression to 
include inhomogeneous strain by expanding the effective Hamiltonian subspace according 
to the procedure described in [pif . The explicit expression is given in PDJ . 



3. 3. Determination of parameters from first principles 

With the truncations described in the previous section, the effective Hamiltonian 
for PbTiOz contains a total of 21 parameters to be determined from first- principles cal- 
culations. These are performed with the preconditioned conjugate-gradients plane-wave 
pseudopotential method, using the program CASTEP 2.1 p3[j, and with the variational 



DFT perturbation method, using a program based on the formalism in |Tl|]. Full details 
of the calculations are provided elsewhere [|2(| . 
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The determination of the quadratic coupling parameters is greatly facilitated by the 
use of the DFT perturbation method for the direct calculation of force constant matrices 
at arbitrary q, circumventing the need for computationally intensive supercell calculations. 
In addition, this method can be used for the direct calculation of Born effective charges Z* 
and the dielectric constant eoo, from which the dipole-dipole term in T~C e ff is immediately 
obtained. At a given q, the energies of the three eigenvectors of the quadratic part of the 
effective Hamiltonian are obtained by computing the quadratic part of the first-principles 
total energy of the ionic displacement patterns in the Aq subspace from the DFT force 
constants at that q. The eight short-range intercell interactions ol, ay, &z,, &ti, cl 
and ct are then obtained by requiring that the expressions for the effective Hamiltonian 
eigenenergies in terms of Z* , and the short range intercell interactions reproduce these 
first-principles values for selected modes at a set of q. This process is illustrated in fig. [| 
The validity of the dipole-dipole form for intercell interactions beyond third neighbor can 
examined by comparing the values predicted by H e f / for calculated energies not included 
in the determination of the parameters. In fig. |3|, these comparisons for q along (111) are 
shown to be very good. 

Given the quadratic coupling parameters, the higher order terms in the onsite potential 
(5, C, D and E), the elastic constants (Cn, C12 and C44) and the couplings between 
strain and local distortions (go, g\ and g^) can be readily determined from the energies 
of configurations in which £3 is independent of i, requiring only total energy calculations 
for primitive unit cells containing five atoms. Further details of this procedure and the 
complete set of numerical parameter values for the PbTiOs effective Hamiltonian are given 
in pO. 



4. Monte Carlo Study of the PbTiOs transition 

The form of the effective Hamiltonian, while greatly simplified relative to the orig- 
inal lattice Hamiltonian, is still sufficiently complicated to discourage the application of 
analytical statistical mechanics methods such as renormalization group analysis or high- 
temperature expansions. However, it is quite suitable for Monte Carlo simulation, since 
changes in energy with configuration are readily calculated numerically |[24| . We use a 
single-flip Metropolis approach, with runs ranging between 25,000 and 150,000 sweeps 
through the lattice. Periodic boundary conditions are imposed on an L x L x L simulation 
cell, with finite-size scaling applied for L ranging from 5 to 11. 



7 



In a finite-size simulation, a first-order transition such as the ferroelectric transition 
in PbTiOz leads to the coexistence of two distinct phases, separated by an energy barrier, 
in a temperature range near T c . Recently developed Monte Carlo methods for first-order 



transitions |?5| determine the transition temperature from the condition that the difference 
in the free energies of the two phases be zero. If the two phases are sampled ergodically 
in the simulation, this difference is quite easy to compute. However, in cases of strongly 
first-order transitions, unbiased sampling of the two phases can be extremely difficult to 
achieve, especially as the simulation-cell size increases. While we are presently considering 
various algorithms suitable for this situation, for now we put bounds on T c by monitoring 
the sensitivity of the average structural parameters to the choice of initial state: T> is 
the lowest temperature at which the system averages are characteristic of the cubic state, 
starting with an inital ground state tetragonal configuration, while T< is the highest tem- 
perature at which a starting cubic configuration results in a tetragonal state. These bounds 
are plotted in fig. |4|. A value of T c = 660-KT, obtained from averaging the bounds at the 
largest system size, is in very good agreement with the experimental transition tempera- 
ture 763K. An estimate of a latent heat of 3400 J/mol, extracted from separate runs for 
the two phases, obtained by using different inital states at T c , is in rough agreement with 
the measured value of 4800 J/mol [[26], and very much larger than the 209 J/mol latent 
heat of the cubic-tetragonal transition in BaTiOs |27] . 



One of the opportunities offered by this first-principles analysis is to investigate the role 
of different contributions in determining the behavior at the transition. One particularly 
striking result of this procedure is that if the strain coupling parameters go, g\ and gi are 
set to zero, the character of the transition changes significantly. Specifically, the system 
is then found to undergo a second- order transition from the high-temperature cubic phase 
to a low-temperature rhombohedral phase with a lower transition temperature of 400K. 
Thus, not only is the strain coupling responsible for stabilizing the ground-state tetragonal 
structure, as has been previously noted 0, but it is also the key factor in producing the 
strong first-order character of the transition. 

Another advantage of computer simulation of a realistic first-principles model is access 
to microscopic information about short-range order. This is particularly of interest for the 
high-temperature paraelectric phase. A microscopically nonpolar character for this phase, 
where the probability distribution of local distortions has a single peak centered at zero, 
is generally associated with a displacive nature for the transition. In general, a single-well 
potential will result in a single-peak distribution if intercell correlations are unimportant, 
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as in the limit of high temperatures. A connection between a single-well onsite potential 
and displacive behavior has been further supported by a dynamical study of a $ 4 model 

At first glance, the situation in PbTiOs appears to be straightforward. The onsite 
potential has a single- well character (positive quadratic coefficient in (|3.5|)), which by 
the above arguments seems to explain the reported displacive character of the PbTiOs 
transition |]J. However, there is recent evidence that this transition cannot be simply 
described as displacive. Recent EXAFS measurements [[28] show that local distortions, 
characteristic of the low-temperature structure, persist at least as high as 85K above T c . 
In our theoretical approach, we can suggest a possible explanation. In the realistic H e ff, 
unlike the $ 4 model, we find that unstable modes of the perovskite structure extend all 
the way to the zone boundary (fig. §). The associated lowering of the energy for local 
distortions, even in the absence of long range order, could significantly affect the short- 
range order above T c . This feature, rather than the sign of the quadratic coefficient in 
the onsite potential, may be the signature of a microscopically polar paraelectric phase. 
More detailed Monte Carlo simulations, including the effective Hamiltonian generalized to 
include inhomogeneous strain, are in progress. In addition, effective-Hamiltonian molecular 
dynamics studies are planned which could reconcile large local distortions with the observed 
"displacive" behavior of the soft mode above T c . 



5. Further applications 

The lattice Wannier function method for the construction of effective Hamiltonians 
is applicable to finite-temperature structural properties in a variety of systems. Effective 
Hamiltonian studies of structural phase transitions currently in progress include the anti- 
ferroelectric transition in PbZrO^ ]2£j , the ferroelectric transition in KNbOs [ |30f , and the 
cubic-tetragonal transition in ZVO2 [31]]. Effective Hamiltonians can also be constructed 
to model the temperature dependence of structural parameters in systems like pyroelectric 
BeO Because of the simplicity of the kinetic energy in a Wannier basis constructed 
from the dynamical matrix eigenvectors, this approach is readily extended beyond the 
equilibrium classical statistical mechanics analysis we have presented here. For example, it 
is straightforward to construct a quantum mechanical model, for example, for a quantum 
paraelectric such as SrTiO^ |33| . In addition, important issues in the dynamical properties 
of ferroelectrics, such as soft-mode behavior and switching, can be studied in molecular 
dynamics simulations for individual materials. 
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6. Summary and Conclusions 

In conclusion, we have described a systematic method for the construction of effec- 
tive lattice Hamiltonians from first-principles calculations through the use of a localized, 
symmetrized basis set of "lattice Wannier functions." Applying this method to the fer- 
roelectric transition in PbTiOs, we find a first order cubic-tetragonal transition at 660 
K. The involvement of the Pb atom in the lattice instability and the coupling of local 
distortions to strain are found to be particularly important in producing the behavior 
characteristic of the PbTiO^ transition. A tentative explanation for the presence of local 
distortions experimentally observed above T c is suggested. Further applications of this 
method to a variety of systems and structures are proposed for first-principles study of 
finite-temperature structural properties in individual materials. 
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Figure Captions 

Fig. 1. A sketch showing the features of the dispersion relation ui\(q) characteristic of 
systems undergoing structural transitions, and a possible decomposition into the 
two subspaces A s and Ao described in the text. 

Fig. 2. Independent displacement patterns for the innermost coordination shells of the 
lattice Wannier functions transforming like the z vector component at the Pb 
atom positions. The open squares represent the Pb ions, the shaded squares the 
Ti ions, and the open and shaded circles represent inequivalent oxygen ions. The 
three additional displacement patterns for first and second neighbor Pb shells, in- 
cluded in the analysis, are not shown. The x- and y-like components are obtained 
by rotating the patterns shown by ^ around an appropriate cartesian axis. 

Fig. 3. Determination of quadratic short-range intercell interaction parameters for 
PbTiO^. The solid and open circles indicate the first-principles values for the 
energies of effective Hamiltonian eigenvectors. The solid circles are used to deter- 
mine values for the short-range interaction parameters, which are used to generate 
the effective Hamiltonian dispersion curves, shown as solid lines. The open cir- 
cles are first-principles results not included in the parameter determination, for 
comparison. 

Fig. 4. Monte Carlo estimate of T c as a function of increasing simulation-cell size L. At 
each L, the vertical line extends from the lower bound T < to the upper bound 
T<, calculated as described in the text. 
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